package mg2d.solver;

import mg2d.Matrix;
import mg2d.Tensor;
import mg2d.domain.Domain;
import mg2d.equation.ConvDiffEqn;

public abstract class Multigrid {
	Matrix[] f;
	Matrix[] u;
	Matrix[] r;
	Tensor[] coeff;
	int n;
	int levels;
	double h;
	ConvDiffEqn cde;
	Domain domain;
	
	public Multigrid(ConvDiffEqn cde, Domain domain, int n) {
		this.cde = cde;
		this.domain = domain;
		this.n = n;
		this.levels = log2(n);
		discret(levels);
	}
	
	/*
	 * Define the functions for the convection diffusion equation
	 * u_xx + u_yy + p*u_x + q*u_y = f
	 */
	
	/**
	 * The analytic function f in u_xx + u_yy + p*u_x + q*u_y = f.
	 * @param x
	 * @param y
	 * @return The value of f at the point (x,y).
	 */
	double f(double x, double y) {
		return cde.f(x, y);
	}
	double f(int i, int j, int level) {
		return f(x(i,j,level), y(i,j,level));
	}
	/**
	 * The analytic function u in u_xx + u_yy + p*u_x + q*u_y = f.
	 * @param x
	 * @param y
	 * @return The value of u at the point (x,y).
	 */
	double u(double x, double y) {
		return cde.u(x, y);
	}
	double u(int i, int j, int level) {
		return u(x(i,j,level), y(i,j,level));
	}
	/**
	 * The analytic function p in u_xx + u_yy + p*u_x + q*u_y = f.
	 * @param x
	 * @param y
	 * @return The value of p at the point (x,y).
	 */
	double p(double x, double y) {
		return cde.p(x, y);
	}
	double p(int i, int j, int level) {
		return p(x(i,j,level), y(i,j,level));
	}
	/**
	 * The analytic function q in u_xx + u_yy + p*u_x + q*u_y = f.
	 * @param x
	 * @param y
	 * @return The value of q at the point (x,y).
	 */
	double q(double x, double y) {
		return cde.q(x, y);
	}
	double q(int i, int j, int level) {
		return q(x(i,j,level), y(i,j,level));
	}
	/**
	 * Convert index i to x coordinate.
	 * @param i
	 * @param j
	 * @param level
	 * @return The x coordinate at index i.
	 */
	double x(int i, int j, int level) {
		return domain.x(i, j, level);
	}
	/**
	 * Convert index j to y coordinate.
	 * @param i
	 * @param j
	 * @param level
	 * @return The y coordinate at index j.
	 */
	double y(int i, int j, int level) {
		return domain.y(i, j, level);
	}
	
	/**
	 * Populate u (boundary), f, and coeff.
	 * @param levels
	 */
	abstract void discret(int levels);
	abstract void restriction(int fineLevel, int coarseLevel);
	abstract void prolongation(int coarseLevel, int fineLevel);
	abstract double relax(int level);
	abstract double coarseRelax();
	abstract double residual2norm();
	
	public double[][] solve(double tolerance) {
		double L2Norm = 1;
		int iter=0;
		double r0 = 1;
		double rn = 1;
		while(L2Norm > tolerance && iter < 10000) {
			for(int l=0; l<levels-1; l++) {
				L2Norm = relax(l);
				//System.out.println(l+" "+L2Norm);
				restriction(l, l+1);
			}
			
			L2Norm = coarseRelax();
			//System.out.println(levels+" "+L2Norm);
			
			for(int l=levels-1; l>0; l--) {
				prolongation(l, l-1);
				L2Norm = relax(l-1);
				//System.out.println((l-1)+" "+L2Norm);
			}
			//printError();
			if(iter==0) {
				r0 = L2Norm;
			} else {
				rn = L2Norm;
				L2Norm = rn/r0;
				//System.out.println(L2Norm);
			}
			iter++;
		}
		System.out.println();
		System.out.println("V-cycles completed.");
		System.out.println(iter+" cycles to converge.");
		//System.out.println("Norm: "+L2Norm);
		
		return u[0].m;
	}
	
	public static int log2(int n){
	    if(n <= 0) throw new IllegalArgumentException();
	    return 31 - Integer.numberOfLeadingZeros(n);
	}
	
	public int levelN(int level){
		if(level==0)
			return n;
		return n/(2 << (level-1) );
	}
	
	public double printError() {
		double maxError = 0;
		double error;
		int maxi = 0;
		int maxj = 0;
		for(int i=0; i<=n; i++) {
			for(int j=0; j<=n; j++) {
				error = Math.abs(u[0].m[i][j] - u(i,j,0));
				if(error > maxError) {
					maxError = error;
					maxi = i;
					maxj = j;
				}
			}
		}
		System.out.println("Maximum error: "+maxError+" at ("+maxi+", "+maxj+")");
		//System.out.println("Value of u at point: "+u(maxi,maxj,0)+" "+u[0].m[maxi][maxj]);
		return maxError;
	}
	
	public void test() {
		System.out.println("f at (0,0): "+f(0,0));
		System.out.println("f at i=0,j=0: "+f(0,0,0));
		System.out.println("f at i=0,j=0: "+f(0,0,2));
		System.out.println();
		System.out.println("u at (0,0): "+u(0,0));
		System.out.println("u at i=0,j=0: "+u(0,0,0));
		System.out.println("u at i=0,j=0: "+u(0,0,2));
		System.out.println();
		System.out.println("p at (0,0): "+p(0,0));
		System.out.println("p at i=0,j=0: "+p(0,0,0));
		System.out.println("p at i=0,j=0: "+p(0,0,2));
		System.out.println();
		System.out.println("q at (0,0): "+q(0,0));
		System.out.println("q at i=0,j=0: "+q(0,0,0));
		System.out.println("q at i=0,j=0: "+q(0,0,2));
		System.out.println();
		System.out.println("u at (mid,mid): "+u(Math.sqrt(3.0)/6.0,.5));
		System.out.println("u at i=32,j=32: "+u(n/2,n/2,0));
		System.out.println("u at 1,1: "+u(1,1,levels-1));
	}
	
	public void testRelax() {
		int iter = 0;
		double initialNorm = relax(0);
		iter++;
		double norm = 1.0;
		while (norm > 1e-10 && iter < 10000) {
			iter++;
			norm = relax(0);
			norm = norm/initialNorm;
			//System.out.println("GS iteration: "+iter+" norm: "+norm);
			//printError();
		}
		System.out.println("iterations: "+iter);
		printError();
	}
}
